Novel use of online optimization in a mathematical model of COVID-19 to guide the relaxation of pandemic mitigation measures

Since early 2020, non-pharmaceutical interventions (NPIs)—implemented at varying levels of severity and based on widely-divergent perspectives of risk tolerance—have been the primary means to control SARS-CoV-2 transmission. This paper aims to identify how risk tolerance and vaccination rates impact the rate at which a population can return to pre-pandemic contact behavior. To this end, we developed a novel mathematical model and we used techniques from feedback control to inform data-driven decision-making. We use this model to identify optimal levels of NPIs across geographical regions in order to guarantee that hospitalizations will not exceed given risk tolerance thresholds. Results are shown for the state of Colorado, United States, and they suggest that: coordination in decision-making across regions is essential to maintain the daily number of hospitalizations below the desired limits; increasing risk tolerance can decrease the number of days required to discontinue NPIs, at the cost of an increased number of deaths; and if vaccination uptake is less than 70%, at most levels of risk tolerance, return to pre-pandemic contact behaviors before the early months of 2022 may newly jeopardize the healthcare system. The sooner we can acquire population-level vaccination of greater than 70%, the sooner we can safely return to pre-pandemic behaviors.

A prototypical convex optimization problem can be expressed in the form u * ∈ arg min u∈R n f (u), (S.3a) subject to: u ∈ U, (S.3b) where f : R n → R is a convex function and U ⊆ R n is a convex set. The set U is called feasible set. Problem (S.3) is feasible if the set U is not empty. We say that u * is a global minimum of f over U if f (u) ≥ f (u * ) for all u ∈ U. An optimization problem in functional form is a problem of the form (S. 3) where U is specified as the intersection of convex sets, namely, U := {g i (u) ≤ 0 and h j (u) = 0, i = 1, . . . , m, j = 1, . . . , p}, where g i : R n → R are convex functions and h i : R n → R are affine functions.  A key result in convex optimization [2,Thm. 9.2] states that if f : D → R is a convex and continuously differentiable function and D ⊆ R n is a closed and convex set, then u * is a global minimum of f over U if and only if u * is a stationary point of (S. 3). Differently, if f is not convex, then stationary points can be local or global minima, maxima, or saddle points.
By building upon the presented convex optimization framework, we next discuss some properties of the feedback controller (4) utilized in this work. Let f : D → R be a continuously differentiable function and let U ⊆ R n be a closed and convex set. In (4), we we use the following dynamics as a feedback controller for (1):ẋ = P U u − η∇f (u) − u, (S. 4) where η > 0 is a design parameter. The solutions of the dynamical system (S.4) satisfy the following properties (proven in [3,4,5]).
ii) If f has Lipschitz-continuous gradient in its domain, then for each u 0 ∈ D there exists a unique continuously differentiable solution u(t) to (S.4) with u(0) = u 0 that is defined for all t ≥ 0.
i) If u 0 / ∈ U, then u(t) approaches the set U exponentially-fast. Furthermore, the set U is forward invariant.
ii) If f is convex, then every stationary point of (S.3) is a Lyapunov stable equilibrium point of (S.4).
iii) If f is strongly convex, then every solution of (S.4) converges to the unique solution of (S.3) exponentially fast.
Finally, we note that if f is not convex, then (S.4) converges to the a stationary point of the problem (S.3); however, only strict local optimizers are asymptotically stable [6].

Supplementary Note 2: Costs and estimated constraints
As discussed in the main paper, the proposed technical approach builds upon formulating an optimization problem that captures the desired social and economic metrics, and incorporates constraints related to the largest admissible number of hospitalizations. In our work, we select cost functions φ i : [0, 1] → R ≥0 of the decision variable u i in region i that model the societal cost induced by the introduction of NPIs in the region, including the economic impact of the raised control measures, and/or the societal response to restriction orders. Mathematically, we require that u → φ i (u) is convex and differentiable function with a Lipschitz-continuous derivative. Since societal losses do not increase as NPIs are lifted, we also assume that u → φ i (u) is a non-increasing function in its domain. Relevant examples of such a functions include: where c > 0 is a free parameter. The Example 1 above was proposed in [7], where the choice of c is made according to the number of employees that are affected by NPIs in the region i. Additional functions are discussed in [7], and additional remarks on the cost of NPIs are provided in e.g., [8,9]. In all the numerical analyses presented in the paper, we used the function Example 2 with c = 1.
To capture the relationship between control variables and the number of infectious individuals at the endemic equilibrium (i.e. when t = ∞), we denote by u → F i (u) the function that maps the instantaneous values of NPIs u across regions to the fraction of infections in the region i at the endemic equilibrium. For the single-region model (3), a closed-form expression for the this map is available (see Supplementary Note 3); however, for the multi-region model (4), the map F i (u) for each region i = 1, . . . , N is not available in closed form. To estimate each of the maps F i , several data-based approaches can be used by using the epidemic transmission model as a predictor. In particular, for all the numerical simulations presented in the paper, we simulated the dynamics (4) with initial condition given by instantaneous model state and for different values of the control parameter {u (1) , . . . u (p) }. Then, we used use the output of the simulations (precisely, the values of the steady-state infections evaluated at at the trial points {u (1) , . . . u (p) }) to construct an approximate map u → F i (u). Alternative functional methods that can be used include Gaussian Processes, Kernel-based Regression, Strongly Convex Regression, or Linear Regression. In the main paper, we approximate the map F i (u) using a linear function of the formF i (u) = a i u + b i , for a i ∈ R N and b ∈ R. In particular, the p test points {u (1) , . . . u (p) } are chosen in a neighborhood of the current vector u, and we run the model (4) to obtain the data {(u (1) , F i (u (p) )), . . . , (ū (p) , F i (u (1) ))}, where the model is run with initial state given by the current epidemic state x(0) = x and with constant levels of NPIs {u (1) , . . . u (p) }. With this data set, define: where the subscript x 0 is used to highlight the dependence on the model state x. Then, a linear approximation of F i (u) can be found as follows: The estimate of the map can be updated to reflect changes in the state is a local linear approximant to the true map.
In the main problem formulation, we utilize a function ı * i → ψ i (ı * i ) that models a cost associated with the infections in region i at the endemic equilibrium; we assume that ψ i is convex and differentiable function with a Lipschitz-continuous derivative. Examples include: where k ≥ 0 is a parameter and i th is a given threshold associated with risk tolerance. The parameter k can be used to model the cost incurred by the healthcare system at a given number of infectious individuals (cost of testing, hospitalizations, etc) or economic costs related to possible quarantines and missing work days. SinceF i (u) is linear, we have the the composition of functions u → ψ i (F i (u)) is convex and differentiable, and it has a Lipschitz-continuous gradient [2] (under the condition that a i < ∞ by design). If methods such as Gaussian Processes are utilized to estimate the map F i , then convexity of F i (and, hence, of ψ i (F i (u))) is not guaranteed. In this case, the convergence and stability properties of gradient flow methods change; in particular, the asymptotic stability of the gradient flow explained in Supplementary Note 1 is a property that holds only for strict local minima. The technical approach explained in the paper also relies on the map u → H i (u; x) for the region i, with which is a function that maps the instantaneous level of NPIs u (in all the regions) into the peak of hospitalizations in the region i, given the current state of the model. A closed form expression for this map is generally not available; we adopt simulation-based approach similar to the one explained above for the maps F i , and we utilize a linear approximation of the formĤ i (u; x) = m i u i + q i (see Supplementary  Fig.2). In this case, the approximate constraint in the problem (3) becomes m i u i + q i ≤ h lim,i , and it is convex (see Supplementary Note 1).
We conclude by noting that in this work we utilize the Susceptible-Exposed-Infectious-Hospitalized-Recovered-Vaccinated-Susceptible (SEIHRVS) transmission model, explained in the section Methods. However, alternative networked models could be utilized for the simulation-based estimation of the maps (for example, a network version of the recently-proposed model SIDARTHE [10]).

Supplementary Note 3: Single-region model
In this section, we provide details on the single-region model and on the implementation of the related feedback controller. We recall that the single-region model reads (see (1)): where we recall that u ∈ [0, 1] specifies the level of permitted social activity. In this case, an exact expression for F can be derived in closed-form, which reads as: .

(S.7)
We notice that, when y = 0 (i.e. the model does not incorporate vaccinations) and σ = 0 (i.e. there is no loss of immunity), the expression (S.7) simplifies to the well-known input-to-steady state map valid for the SEIR model [11,Section 9.3]: where R 0 := βu ( +δ)(γ+δ) denotes the basic reproduction number. For this single-region model, the problem of optimizing the severity of NPIs while maintain the number of COVID-19 hospitalizations under a pre-specified threshold (denoted by h lim ) can be formulated as follows: (loss depends on lockdown policy and infections at steady-state) where u → H(u; x) is a function that maps NPIs into peak hospitalizations. Since F(u) is increasing in u and ψ(ı) is convex in its argument, we have that ψ(F(u)) is convex.
The function H(u; x) can be approximated numerically. In particular, a linear approximationĤ(u; x) can be found as described in Supplementary Note 3. Alternatively, we also point out that numerical evidences suggest that the map H(u; x) is monotonically increasing in the variable u, for a fixed x; in this case, the constraint H(u; x) ≤ h lim can be converted into u ≤ u lim , where u lim is such that H(u lim ; x) = h lim .
Define the setÛ x := {u : u ∈ [0, 1] andĤ(u; x) − h lim ≤ 0} (alternatively,Û x := [0, min{1, u lim }] if H(u; x) is monotonically increasing in the variable u). The controller proposed in the main paper, in this case, reads as follows:u where η > 0 is a tunable parameter of the controller. We note that the controller (S.9) leverages two types of feedback from the system: (i) it uses the instantaneous fraction of infectious individuals ı; (ii) it relies on a projection onto the setÛ x , which is parametrized by the instantaneous state of the system.
For these reasons, the control dynamics (S.9) describe a dynamic state-feedback controller for the epidemic model (S.6). A block diagram of the controller is illustrated in Supplementary Fig. 1.
We conclude by discussing a special case in which the number of infections at the endemic equilibrium is disregarded from the loss function. In this case, (S.8) can be simplified to optimize the severity of the NPIs while maintaining the hospitalizations under a given limit as follows: (cost depends on lockdown policy) Moreover, we note that, when φ(·) is a monotonically decreasing function in its domain, an optimal NPI policy that solves (S.10) can also be computed in closed form by letting: That is, the restrictions on social interactions are lifted to a level where the forecasted amount of hospitalizations will be precisely at the limit h lim . Although in this case the solution u * can be found explicitly, one can still utilize the controller in (S.9) to adopt a less aggressive approach for lifting the NPIs. In particular, the controller in (S.9) is simplified in this case aṡ where the setÛ x is re-estimated on a daily basis based on the current number of infections, vaccinations, and hospitalizations.

Supplementary Note 4: Stability and convergence of the control method
We next discuss the stability and convergence of the NPI controller when coupled with the dynamics of epidemic transmission.
Projected gradient flow. Based on the discussion in Supplementary Note 1, when the set U is constant over time and the cost function u → φ(u) + ψ(F(u)) is convex, continuously differentiable, and has a Lipschitzcontinuous gradient over the domain of the function, the projected gradient floẇ where J(u) is the Jacobian matrix of F evaluated at u, is globally asymptotically stable [3] (in the sense that every trajectory u converges to the set of minimizers); on the other hand, (S.12) is globally exponentially stable when the cost is strongly convex over the set U [3], in the sense that every trajectory u converges to the unique minimizer u * exponentially fast (see also [4,5] for a proof of these claims) In our case, the controller relies on the approximate setÛ x ; this set is a polyhedron (it is given by the intersection of finitely many half-spaces) and it is convex. The polyhedronÛ x continuously evolves over time since it is parametrized by the state x(t) of the SEIHRVS model. Let u * denote an optimal solution of the optimization problem min and assume that u * is unique for simplicity of exposition. Under the assumption that the map t →Û x is locally Lipschitz (and therefore the projection onto moving polyhedra is locally Lipschitz [12]), the solution trajectory u(t) ofu = PÛ x u − η(∇φ(u) + J(u) ∇ψ(F(u))) − u exists, it is unique, and it can be shown to converge to the optimizer of (S.8) (see [3]).
Dynamical transmission model. Consider the single-region Susceptible-Exposed-Infectious-Hospitalized-Recovered-Vaccinated (SEIHRV) model, which corresponds to the dynamical system (2) in the main paper when the loss of immunity is discarded from the model. The SEIHRV dynamical model has an endemic equilibrium that is globally asymptotically stable for any u; this can be shown by extending the arguments in [11,Section 9.3] for the single-population Susceptible-Exposed-Infectious-Recovered (SEIR) model. The networked SEIHRV model has a unique endemic equilibrium. Assuming that the reproduction number is strictly larger than 1 and the interaction graph (modeling mobility) is strongly connected, it can be shown that the endemic equilibrium is globally asymptotically stable [11,Section 9.3].
We conclude this part by noting that the stability of the full dynamical model with the loss of immunity is an open problem; though, we also point out that it is still unclear how to model the loss of immunity of COVID-19 vaccines (the simulations in this work are run for the period of one year, under the assumption that there is no loss of immunity in this time frame).
Interconnection between controller and dynamics of the epidemic. The stability of the interconnection between a gradient-flow controller and a dynamical system can be investigated using arguments from singular perturbation theory or small-gain stability [13].
Consider the single-region SEIHRV dynamical model (1) interconnected with the projected gradient floẇ In this particular case, we have an interconnection between the SEIHRV dynamics, which have an endemic equilibrium that is globally asymptotically stable, and the projected gradient flow. In this setting, it is possible to utilize the procedure explained in [13,14] to show that there exists a rage of values for the gain η of the projected gradient flow that renders the optimizers of the optimization problems input-to-state stable. Similar arguments can be utilized for the interconnection of the multi-region SEIHRV dynamics with the projected gradient flow controller in equation (1) of the main paper.
We conclude by noting that an open problem is to establish stability of the multi-region SEIHRV dynamical model interconnected with the distributed controllers in equation (2) in the main paper; see also Supplementary Fig. 1. The analysis of this interconnection calls for new theoretical results that are not covered by the works of e.g., [3,6,13,14].

Supplementary Note 5: Data
We organized the model-fitting phase into two main stages. First, we fitted the SEIHRVS model by combining model parameters from [15] with hospitalization data by using a via a prediction-correction algorithm (see table in the main paper). Second, we used cell-phone usage from SafeGraph (https://docs.safegraph. com/docs/social-distancing-metrics) to estimate the interaction matrix of the network.
Precisely, model parameters were selected (via a prediction-correction algorithm) to fit hospitalization data for the period from 1/18.21 to 2/28/21. State-wide hospitalization data was obtained from EMResources and regional-level data came from Colorado COVID Patient Hospitalization Surveillance [16,15]. This timeinterval was selected due to a steady decline in hospitalizations observed over this period, thus overcoming the need to account for varying transmission rates during model-fitting. Validation was carried out by comparing the evolution of the parametrized model with hospitalization data, and by assessing performance through a mean squared error. For the parametrized model, the mean squared error associated with hospitalizations and infections was less than 10% over the entire dataset. The model was initialized on 1/18/21 by using infection data from [17]. The state of the pandemic on 03/01/21 was computed by initializing the model on 1/18/21 and by simulating the evolution of the epidemic with constant vaccination rate of y = 12, 000 vax/day and vaccination efficacy ν = 0.9.
Cell-phone data usage was used to estimate the mobility matrix. Data was obtained from SafeGraph, in the form of counts of cell phone accounts that normally reside in county i and traveled to county j on a specific date, have been averaged over the time interval 03/01/20-12/31/2020.
All computational analyses and the fitting of data were performed using MATLAB and the corresponding optimization toolbox. Our Montecarlo simulations suggest that the outcomes of our simulations are robust to variations of the parameter, confirming their validity to describe reality. Note that the model describing the epidemic spread is highly nonlinear and therefore potentially sensitive to parameter perturbations.  1 Q y I z 8 R V u r j c r P W d G r F b T h i r f 4 F s L U a L R 6 T E F 2 L u U 5 X 0 i X R F 1 f L B O r t M p P P y K u V B l E g W W M u J n M R r y b C V 3 3 v L 5 j G z p D d X D W L F X F 1 r y 3 J J T C y p n l B p F h Y J h 0 / a + a / 2 4 v r a M G / b 4 p K p c W 6 1 W f 5 / 6 V b T u 5 h E 6 r j y I P X L / S T 2 y g P 7 v V 5 P x i S Q Y e i J M j E / c t l t e W w m V I u 1 4 S G 1 4 V J F Q t W v J G Z t M V d n C r x 5 V r 0 j G U a 2 U x 6 V X L 2 Y 6 t D N f d Z o N M Y R i X l g q z f c 6 r C Z 6 h P b l q H H g o l 0 0 7 E 6 l 0 / i C Q + y 9 I f O R v e F a j T G N n P G e j E K G T N p u W l n 4 + e X W f V g N p M u 4 x N X Z q k + u D H 2 m C O X 5 + R q z b 1 g f m N 5 R i v 0 8 3 W 2 P G z x M t K Y 2 R n M p 0 5 c F O o l L M v P F r C 7 / E A S K E 7 k H c 8 u u 1 f p 2 K f h L B 3 T 0 L P z l 9 l a X + u u 5 y c q h n t Z q n J K K U 1 7 W V W t 4 T T 3 f M L Y V w t 6 i h J T E h c S I 9 W r J s w 4 p I W I i R N q N V i s E B n g i F o v N 4 X E o l v N 2 E w t 7 E J I 9 X G G k 0 k x k n e r G V o J U Z 2 q P J u d P G M R r 7 W D H t u B p g N E f a C 0 j 5 / 2 W 3 3 Y W 0 S m J h O R 0 C Q Q 0 Y X l C + I h U Y 3 s 6 q N 3 E W 1 q 2 k R 0 o u k E 0 b G m Y 0 Q D T Q N E 2 5 q 2 E f U 0 9 R D t a 9 p H 9 F r T a 0 R b m r b w e t c 0 R H S o 6 R D R n q Y 9 R C N N I 0 R n m s 4 Q n W s 6 R / R O 0 z t E F 5 o u E J 1 q O s V r x C y v E X O x 1 C t b w s o z 1 G l Z a L M Q E I K E g l C 8 w U B s J A y E 4 b M 5 Q A 6 i C c g E i Q v i I g l B Q i Q R S I T k F u Q W S Q w S I x E g a G / S G c g M S Q K S I L k D u U M y B Z k i 4 S A c y T 3 I P Z I 5 y B z J J g g q C X Q L B G 0 q u g 2 C t j b d A U E V l P Z A 0 J 6 n u y C o X t E D E F R 0 6 R 4 I 2 p / 0 D c g b J H 2 Q P p J D E F Q I l l t J C S r S 9 B g E F U 1 Q y I z 8 R V u r j c r P W d G r F b T h i r f 4 F s L U a L R 6 T E F 2 L u U 5 X 0 i X R F 1 f L B O r t M p P P y K u V B l E g W W M u J n M R r y b C V 3 3 v L 5 j G z p D d X D W L F X F 1 r y 3 J J T C y p n l B p F h Y J h 0 / a + a / 2 4 v r a M G / b 4 p K p c W 6 1 W f 5 / 6 V b T u 5 h E 6 r j y I P X L / S T 2 y g P 7 v V 5 P x i S Q Y e i J M j E / c t l t e W w m V I u 1 4 S G 1 4 V J F Q t W v J G Z t M V d n C r x 5 V r 0 j G U a 2 U x 6 V X L 2 Y 6 t D N f d Z o N M Y R i X l g q z f c 6 r C Z 6 h P b l q H H g o l 0 0 7 E 6 l 0 / i C Q + y 9 I f O R v e F a j T G N n P G e j E K G T N p u W l n 4 + e X W f V g N p M u 4 x N X Z q k + u D H 2 m C O X 5 + R q z b 1 g f m N 5 R i v 0 8 3 W 2 P G z x M t K Y 2 R n M p 0 5 c F O o l L M v P F r C 7 / E A S K E 7 k H c 8 u u 1 f p 2 K f h L B 3 T 0 L P z l 9 l a X + u u 5 y c q h n t Z q n J K K U 1 7 W V W t 4 T T 3 f M L Y V w t 6 i h J T E h c S I 9 W r J s w 4 p I W I i R N q N V i s E B n g i F o v N 4 X E o l v N 2 E w t 7 E J I 9 X G G k 0 k x k n e r G V o J U Z 2 q P J u d P G M R r 7 W D H t u B p g N E f a C 0 j 5 / 2 W 3 3 Y W 0 S m J h O R 0 C Q Q 0 Y X l C + I h U Y 3 s 6 q N 3 E W 1 q 2 k R 0 o u k E 0 b G m Y 0 Q D T Q N E 2 5 q 2 E f U 0 9 R D t a 9 p H 9 F r T a 0 R b m r b w e t c 0 R H S o 6 R D R n q Y 9 R C N N I 0 R n m s 4 Q n W s 6 R / R O 0 z t E F 5 o u E J 1 q O s V r x C y v E X O x 1 C t b w s o z 1 G l Z a L M Q E I K E g l C 8 w U B s J A y E 4 b M 5 Q A 6 i C c g E i Q v i I g l B Q i Q R S I T k F u Q W S Q w S I x E g a G / S G c g M S Q K S I L k D u U M y B Z k i 4 S A c y T 3 I P Z I 5 y B z J J g g q C X Q L B G 0 q u g 2 C t j b d A U E V l P Z A 0 J 6 n u y C o X t E D E F R 0 6 R 4 I 2 p / 0 D c g b J H 2 Q P p J D E F Q I l l t J C S r S 9 B g E F U in a neighborhood of the current epidemic state.